# bcp1.R
# A small example, applying the Bayesian change point functions.
# Bruce Western, April 27, 2004.

source("bcp.R")

x <- c(1,2,3,4,3,5,4,6,5,4,3,7,6,5,8,9,7,0,9,8,7,0,0,3,16,14,13,2,10,13)
x1 <- cbind(1, x)
cp <- 21
n  <- length(x)
ind <- as.numeric(1:n >= cp)
y <- .5 + .2*x  - .3*ind + .2*ind*x + rnorm(n,sd=.5)

pmu <- rep(0,4)
pvb <- diag(rep(100,4))
pa <- .001
pb <- .001

# Analysis with Gibbs sampler
out1 <- gibbs(x1, y, burn=100, iter=2000)

# Trace plots of theta and intercept
plot(out1$theta, type="l")
plot(out1$beta[,1], type="l")

# Plotting posterior for theta
hist(out1$theta)

# Printing p(theta|y), (2000 is number of Gibbs iterations)
cat("Posterior probability of theta (Gibbs)\n")
print(table(out1$theta)/2000)

# Printing posterior expectations and standard deviations for coefs
cat("Posterior means and s.d. of coefficients (Gibbs)\n")
print(apply(out1$beta, 2, mean))
print(sqrt(apply(out1$beta, 2, var)))

# Analysis with SSCP
out2 <- sscp(x1, y, pmu, pvb, pa, pb, bayes=T)

# Plotting posterior for theta
plot(as.numeric(out2[[1]]), out2[[2]][,1], 
     xlab="theta", type="l")

# Printing posterior and log likelihoods for theta
cat("Posterior probability of theta (Chin and Broemling 1980)\n")
print(out2[[2]])

# Printing posterior expectations and standard deviations for coefs
cat("Posterior means and s.d. of coefficients (SSCP)\n")
print(apply(out2$beta, 2, mean))
print(sqrt(apply(out2$beta, 2, var)))



